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Abstract 

Simple approximations for fluxes of atmospheric muons and muon neutrinos are 
developed which display explicitly how the fluxes depend on primary cosmic ray energy 
and on features of pion production. For energies of approximately 10 GeV and above 
the results are sufficiently accurate to calculate response functions and to use for 
estimates of systematic uncertainties. 

1 Introduction 

There are several reasons for which it is useful to have simple approximations for pro- 
duction of muons and neutrinos by interactions of cosmic-ray nucleons in the atmosphere. 
First, the dependence on pion and kaon production, as well as the dependence on energy 
of the primary cosmic rays can be displayed explicitly. This is useful to make estimates 
of the systematic errors in the expected signals in muon and neutrino detectors that arise 
from uncertainties in the primary cosmic-ray spectrum and in meson production, which 
in general depend on energy. 

The dependence of the signal in a particular detector on the primary energy spectrum 
is its response function. The classical use of response functions for muon telescopes |ffl] 
is to interpret measurements of solar modulation of muons as a function of muon energy 
and geomagnetic latitude ||]. In addition, response functions are useful for planning and 
optimizing Monte Carlo simulations of detectors because they display explicitly the relative 
portion of the signal that comes from each part of the primary spectrum. 

For power-law spectra the analytic approximations for fluxes of muons and neutrinos 
are well-known [||, ||, HJ. In this paper we develop approximations which also display the 
dependence on primary energy. Such approximations allow semi-analytic evaluation of 
response functions. They can also be used to estimate the muon-induced signal produced 
by a fluence of protons with an arbitrary energy spectrum and time profile, such as might 
be associated with a strong solar event. 
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The approximate results here are valid for muons and muon-neutrinos at relatively 
high energy, of order ten GeV and above. This includes the interesting cases of under- 
ground measurements of atmospheric muons and of upward energetic muons induced by 
atmospheric neutrinos. The same results can, of course, be obtained with greater precision 
and for low energies, including electron neutrinos, by Monte Carlo simulations or by more 
complicated numerical calculations. The analytic approximations are nevertheless useful 
for the insight they provide and for the diagnostic purposes mentioned above. 



2 Derivation of analytic approximations 

The differential production spectra of muons and neutrinos of energy E from decay of 
pions in atmospheric cosmic ray cascades are given by 

dN S{E,E^) 1 f x dX' 



dE dE n dE N dX ( 



1 - TVjiV d w Jo Atv 

*^rF Nn (E^E N )e- x 'l hN <f>(E N ), 



where P 7T (E 7T , X, X') is the probability that a charged pion produced at slant depth 
X' (g/cm 2 ) survives to depth X (> X'). For muons the kinematic factor, S(E,E n ), 
has the form 

S„(E^o, E v ) = &(E„ - E^o) G(^ - E v ), (2) 
while for muon neutrinos it is 

S^(E u ,E n ) = Q{E V - jt^-t), (3) 

2 

where r w = Here is the Heaviside step function, and the subscript o in Eq. || 

indicates muon energy at production (see Appendix). 

The contribution from kaon-decay has the same form as Eq. |l] with appropriate changes 
in kinematic factors and an overall factor for the branching ratio of Bk-^^u = 0.635. 
The various parameters relevant for N — > [K, w] — * [iv^ are summarized in Table 1. The 
differences between the decay parameters for kaons and pions and between the mass ratios, 
r n and are responsible for the striking difference between the importance of kaons for 
neutrinos as compared to muons. In particular, kaons become the dominant source of 
atmospheric neutrinos for E v larger than ~ 100 GeV even though more pions than kaons 
are produced. 

The integral in Eq. || is over atmospheric depth X' measured in g/cm 2 along the 
trajectory of the incident particle at zenith angle 9. The nucleon interaction length is Ajy. 
The factor exp{— X'/An} (p{E^) is the flux of nucleons of energy En at slant depth X', 
where A^v is the nucleon attenuation length in the atmosphere and 4>(En) is the primary 
spectrum of nucleons per GeV/nucleon. The superposition approximation 0] is used for 
nucleons bound in nuclei. 

The normalized inclusive cross section for iV + air — > -k^ + X is 

F Nw (E„E N ) = d(7( ^' En) * c + (l - xY+ + c_(l - xf-, (4) 
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Table 1: Parameters for atmospheric n + + /i and v„ + 



Mass-square ratios: 








0.573 


.046 


Characteristic Ed eC ay'- 




(-K (-{i 




115 GeV 


850 GeV 1 GeV 


Attenuation lengths: 




Ak Atv Aat 




160 g/cm 2 


180 g/cm 2 120 g/cm 2 86 g/cm 2 


Table 2: Parameters for kaon and pion production 




c+ P+ 


Z + c_ p_ Z_ 


charged pions 


0.92 4.1 


0.047 0.81 4.8 0.034 


charged kaons 


0.037 0.87 


0.0089 0.045 3.5 0.0028 



where x = x T = E^/Ej^. The parameters of this scaling approximation are fixed by 
comparison to accelerator data in the sub-TeV range Q , and the adequacy of the simple 
approximation is verified later by comparison of the calculated neutrino flux to the result 
of a complete Monte Carlo calculation. A more complicated or precise expression could 
be used for the inclusive cross section if desired. The expressions for kaons have the same 
form, and the numerical values used in this paper are given in Table 2. 

The survival probability in Eq. [I] is a product of survival against decay and survival 
against interaction. 

P w (E n ,X,X>) = i^—j exp{-(X-X')/A4. (5) 

Use of the pion attenuation length, A n , rather than the pion interaction length, X n < A n , 
in the interaction factor of Eq. || accounts approximately for the regeneration of charged 
pions when pions interact in the atmosphere. The first factor in Eq. [| is the probability 
for survival against decay, and the explicit form used here is valid in the approximation 
of an exponential atmosphere with scale height ho and decay constant e w = m^hoj ct v . 

Finally, l/d n = t^jiE^X cos 9) is the differential pion decay probability (per g/cm 2 ) 
in an exponential atmosphere. The simple angular dependence displayed here is valid 
when the curvature of the Earth can be neglected {6 < 60°). For large angle the analysis 
is more complicated, but one approach is to use the same form with cos 6 replaced by an 
effective value that depends on 9 0]. 

The next step is to carry out the integral over the production depth, X', in Eq. |j. 
This is done by expanding the exponential of X' and expressing the result as a series of 
the form 

-L / P^,X,X')e- x '^^= C —— £ (-a n )) 1 

a-* Jo Aat Aat ^ VA^ J n\(l + nE n cos9/e 7T ) 

(6) 

where a w = A^/A^v — 1 ~ 0.33. Substituting this expression into Eq. [l] we can now 
integrate over all depths X of pion decay. The upper limit of this integral can safely be 
extended to infinity since most production of secondaries occurs well above the ground. 
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The result for + f M is 



fLA ~ • S " (£ "-^ ) 1 F N „(E w ,E N )G w (E w ,0)<f>(E N ). (7) 



dE v dE^dEjq (1 — r n )E n E n 
The physics of cascade propagation is included in Eq. |7] by the factor 

GAE^e) = ±l £ ij~^ )n '\, ) ■ (8) 

A n , 1 + n E„ cos # /e-r 



n=l 



The contribution from kaons has the same form with ax = Aj^/Ajv — 1 ~ 0.5. 

The dependence on the inclusive cross sections for pion (and kaon) production by nucle- 
ons can be displayed immediately in the approximation of a power-law primary spectrum 
by changing variables from {E u , E-^jEn} to {E u ,x n , z = x^En} and integrating over z to 
obtain 

dN 7 _! F N7r (x n ) f°° dz 

{Xn)' — / —G n {z,9)<j)(z), (9) 



dE v dx^ " 1-rv JE v /(\-r*) 

where 7 is the integral spectral index. The dependence on pion production by nucleons 
appears here as an explicit proportionality factor, the integral of which is the Z-factor, 

Ztv-^tt = / dx (x) 7 " 1 Fjv-*7r(x), etc.. (10) 



The dependence on nucleon production by nucleons and on pion [kaon] production by pions 
[kaons] does not appear explicitly in Eq. |9[ Dependence on these processes is contained 
in the attenuation lengths, An and A n [Ar-]. which appear respectively in Eq. [I] and in 
the expression for G w [Gk] as defined in Eq. ||. The attenuation lengths are related to the 
respective interaction lengths and inclusive cross sections by Aj = Aj/(1 — Z^i), where 
i stands for nucleon (N = p or n), charged pion or charged kaon. Note from Eq. |8| that A,,- 
cancels out of the expression when E^/e^ <C 1, reflecting the fact that low energy pions 
[kaons] usually decay before they interact. The small cross processes <-> 7r ± , as well 
as neutral kaons, are neglected as sources of v and /i in these approximations. 



3 Response Functions and Fluxes 



To obtain the response function, the integral over pion energy is carried out in Eq. |7] 
making use of the scaling approximation of Eq. |j. The flux of Vu + with zenith angle 
6 from decay of tt + + ir~ is given by 



En 



dN„ 



dE u dE]y 



4>(E N ) 



1 



dx G 1T (EpfX, 8) 



(11) 



where x v = E v /En- The contribution from kaons has the same form. 

In the absence of muon decay and energy loss, the expression for the response function 
for muons would be identical to that for neutrinos except that the limits of the integration 
in Eq. 11 are changed to those appropriate for muons (Eq. instead of H). They are 



E t*,o 
En 



and x r 



mm 



1, 



4 




10° 10 1 10 2 10 3 10 4 10 1 10 2 10 3 10 4 
E v (GeV) E v (GeV) 

Figure 1: Comparison of the analytic approximation (solid) for the flux of Vu + with 
the result of the full Monte Carlo (dashed) §§. 



In the Appendix an approximate treatment is given for the effects of muon decay and 
energy loss that can be used for ~ 10 GeV and higher. In this approximation, 



E. 



2GeV\ 
k« + -J 



cos 6 



and 



with 



E N 
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dEucLE^ 
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1 — tv 
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dx G n ^(ENX,9) 



F n (x) 
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E 
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1 + tiEn x j cos 9\le m 



(12) 



(13) 



The functions A{E fjL ) and C n (E fl ) are defined in the Appendix. 

The differential fluxes are then obtained in both cases by summing the contributions 
from pions and kaons and integrating over En. The analytic approximations for the fluxes 
of v^ + Vp and are compared with the results of the full Monte Carlo calculation || 

at two angles in Figs |l| and The values of parameters that have been used in the analytic 
approximations are listed in Tables 1 and 2. They correspond most closely to the hadropro- 
duction model of Ref. |8| for interaction energies around 100 GeV. The approximation falls 
below the neutrino flux at low energy primarily because the contribution from muon decay 
is not included in the approximation. At high energy the approximation is 10-15% higher 
than the Monte Carlo, because the simple, energy-independent forms for kaon production 
used here are not identical to the actual representation in the Monte Carlo. The analytic 
approximations could be tuned to give a better representation of this || or any other 
interaction model, but doing so is not the object of this paper. The approximate results 
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Figure 2: Comparison of the analytic approximation (solid) for the flux of fi + + fj, with 
the result of the full Monte Carlo (dashed) @. 



are sufficiently accurate that the analytic calculation may be used for diagnostic purposes 
above 10 GeV where the contributions from muon decay are unimportant. These include 
propagation of errors, calculation of response functions and evaluation of the phase space 
regions that are most important for the determining the fluxes of neutrinos and muons. 

In Fig. H I show the fractions of + from and from K to demonstrate the 
importance of kaon production. Note that over the entire energy region, approximately 70- 
75% of (f> v comes from free protons (most of the rest is from nucleons in primary helium) Jy]. 
Thus, any uncertainty from the approximate, superposition model treatment Q of bound 
nucleons will be small. 

Having established that the method is adequate, the next step is to calculate the 
response functions for a range of energies. This is illustrated in Fig. ||| for muons and 
neutrinos. Note that the neutrino response is shifted to higher primary energy than the 
response for muons. In addition, it is slower to reach an asymptotic form. These features 
are consequences of the kinematics of it — > fj,Un as compared to K — > [iv^ coupled with the 
greater relative importance of kaons for neutrinos at high energy. Pions must have more 
than twice the energy of daughter u^, while E w > E^q (compare Eqs.^j3|). This shifts the 
Vu response to higher energy when the pion chain is dominant. At higher energy, where 
kaons dominate the production of v„ the response shifts to lower energy. 

4 Estimates of Uncertainties 

To illustrate how the formalism presented here can be used to estimate uncertainties in 
the calculated result, we consider the example of vertically upward, neutrino-induced, 
through-going muons at Super-Kamiokande |lQ| . The distribution of neutrino energies 
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Figure 3: Plot showing fractional contribution of pions and kaons to the flux of v^ + v^. 

that produce this signal (in the absence of oscillations) is shown as the light solid line [|ll]] 
in Fig. ||. The heavy line shows the corresponding distribution of parent nucleon energies. 
This distribution should be used to weight the uncertainties in the input in a convolution 
over primary energy. A similar analysis can be made for neutrino-induced muons that 
stop inside the detector. 

There are two main sources of uncertainties in calculations of fluxes of muons and 
neutrinos with energies in the multi-GeV range and above. We deal first with those 
arising from uncertainties in the primary spectrum and then with those from uncertainties 
in treatment of pion and kaon production. 

4.1 Primary spectrum 

Since approximately 75% of high energy atmospheric neutrinos are produced by incident 
protons, I concentrate on that component of the cosmic-ray spectrum for the analysis of 
uncertainties. A summary of data is shown in Fig. |6| The solid line is the proton spectrum 
used in Ref. ||], which I take as the reference spectrum. The resulting uncertainty in the 
flux of atmospheric neutrinos that is relevant for the Super-K example will be a convolution 
of the uncertainties in the data of Fig. |6| with the response function shown as the heavy 
solid line in Fig. [j| 

It is natural to divide the discussion into three energy regions: 

4.1.1 Up to ~ 100 GeV 

Primaries with energies below 100 GeV are the most precisely measured. They make 
the major contribution to contained neutrino events, but only a small contribution to 
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Figure 4: The distributions of primary energy per nucleon (En) that give rise to vertical 
^u+^u of energy E v (right-hand set of curves) and to fi + + pT . Each curve is normalized to 
unity. For each set the lowest energy is the rightmost curve. The response is approximately 
energy-independent for E v > 3 TeV and for En > 30 GeV when expressed in terms of 
E/En- 



neutrino-induced upward muons. While they do not cover the full range of interest they 
can be used to renormalize the single measurement that spans the important energy 
range up to ~ 1 TeV. Several measurements have been made in the past decade with 
balloon-borne magnetic spectrometers capable of making relatively precise momentum 
measurements for well-identified protons and nuclei. Since the LEAP experiment []12[|, 
five other experiments [13, 14, |l5|, 16, 17] confirm a significantly lower normalization than 
the earlier standard reference |18| . Particularly noteworthy are the BESS98 |l7j and the 
AMS measurements, which agree with each other to within 5% over two decades in 
energy. These data have been used as the basis of a new fit || as shown in Fig. || The 
systematic uncertainty in the primary spectrum below 100 GeV is estimated in Ref. |J as 
< ±5%. 



4.1.2 Up to ~ 1 TeV 

The energy region from 0.1 to 2 TeV, explored by the large balloon-borne tracking calorime- 



ter of Ryan et al. [19], is of great importance for z/^-induced upward muons. We renor- 
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Response for vertically upward muons in Super-K 
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Figure 5: Distribution of neutrino energy (light solid) and primary nucleon energy (heavy 
solid) that produce the through-going, vertically upward signal at Super-Kamiokande. 
Each curve is normalized to unity. The contributions of pions (double dash) and of kaons 
(dotted) are shown as a function of neutrino energy. 



malize the Ryan et al. data downward by 25 % to agree with the magnetic spectrometer 
data below 100 GeV ||. Note that the normalization can reflect a smaller shift in energy 
calibration because of the steep spectrum H. This renormalization brings the measure- 
ment of Ryan et al. into agreement with the reference spectrum. The remaining overall 
uncertainty in this energy region is ~ ±10%. 

4.1.3 The multi-TeV region 

Most of the data in this energy range comes from balloon flights of emulsion chambers, |2Q , 



21, |23T] which measure only the electromagnetic component of the energy deposited by the 
interacting proton in the chamber. This is subject to large fluctuations superimposed on 
the uncertainty in the underlying inelasticity distribution of protons when they interact 
in the emulsion chamber. There is also one measurement in this energy region made with 
a tracking calorimeter carried on flights of the SOKOL satellite [^2] . 

The uncertainties in these measurements are at the level of ±25%. In addition, there 
is an indication that the flux falls below the reference spectrum above ~ 20 TeV by up to 
a factor of 2. However, the contribution from protons with > 20 TeV to the energy 
range responsible for z^-induced upward muons is small. 
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Figure 6: Summary of measurements of the spectrum of cosmic-ray protons. Data are 



from [12, |3L 14, 15, H, 17, |8, 19, 20, 21, 22, 



4.2 Dependence on pion and kaon production 

The dependence on the phase space for pion and kaon production is displayed in Eq. ||, 
which is illustrated for several values of neutrino energy in Fig. |. The fi gure illustrates 
the growing relative importance of kaons. It also illustrates the importance of the extreme 
forward fragmentation region for the contribution of kaons, reflecting the leading p — > 
A + K + process. The shapes are scale invariant in the approximations used here, and 
they indicate the regions of phase space that are important for production of atmospheric 
neutrinos. 

A global estimate of the uncertainty due to pion and kaon production is obtained, as in 
Ref. H, by estimating the uncertainty in the Z-factors for pions and kaons and weighting 
them according to Figs. |3] and [|. 
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Figure 7: Phase space distributions for parents of + from decay of 7r ± and K ± at four 
values of E v . Solid lines are for pions and dashed for kaons. The curves are normalized to 
unity at each energy for the sum of ir — > v and K — > v. from pions and from kaons. 



4.3 Illustration 

To obtain a simple estimate of the overall uncertainty in the signal of upward-moving 
muons due to the uncertainty in the spectrum of atmospheric neutrinos, we combine 
the estimates of several sources of uncertainty as if they were independent, statistical 
uncertainties. For upward, through-going muons in Super-Kamiokande, as shown in Fig. ||, 
the three energy regions of the primary spectrum, < 0.1 TeV, 0.1 < En < 1 TeV and 
> 1 TeV contribute respectively fractions of 0.2, 0.4 and 0.4 to the signal. Weighting and 
adding the corresponding three uncertainties (±.05, ±.10, ±.25) in quadrature gives an 
overall uncertainty of ±11% due to the primary spectrum. Uncertainties in Zjy n and Znk 
are estimated in Ref. [||] as ±12% and ±17% respectively. Roughly equal fractions of the 
signal are from pions and kaons, so the combined uncertainty from the interaction model 
is ~ ±10% giving an overall estimated uncertainty of ±15%. 

For upward, stopping muons in Super-K, the corresponding analysis shows that essen- 
tially all the signal comes from primaries with E < 1 TeV/nucleon, of which ~ 60% are 
in the well-measured region between 10 and 100 GeV. Thus the weighted uncertainty due 
to uncertainty in the primary spectrum is estimated as ±5%. For this lower energy sam- 
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pie, approximately two-thirds are from decay of pions, so the weighted uncertainty from 
hadronic interactions is also somewhat smaller than for through-going muons, ±10%, giv- 
ing a combined uncertainty of ±11%. 

Combining these two numbers results in an estimate of ±19% for the uncertainty in the 
calculated ratio of stopping to through-going neutrino-induced, upward muons in Super-K 
expected in the absence of oscillations. 



5 Summary 

The main results of the paper are contained in Eq. || and in Eqs. 11 and 12 [and the 
corresponding expressions for kaons]. Eq. ||] displays explicitly how the flux depends on 
the phase space for pion [kaon] production. One can see from the corresponding plots 
which regions of phase space need to be measured in experiments such as HARP [gj], in 
order to improve the input to calculations of atmospheric neutrinos. 



Eqs. 11 and 12 show respectively the response functions for Vn and for muons. These 
can be convolved with the response of an underground detector to neutrinos or muons to 
show how the signal depends on the primary spectrum. The response functions in Fig. || 
can also be used as an aid in planning a Monte Carlo calculation of atmospheric muon or 
neutrino fluxes and their signals in a particular detector. By knowing the distribution of 
primary energies that contribute to the neutrino (or muon) flux at , one can determine 
in advance the range and distribution of primary energies needed to obtain a result with 
sufficient statistical accuracy. 

Since the formulas |ll] and 12 display the dependence on primary energy explicitly, 



they can also be used to estimate the signal from a primary spectrum of arbitrary shape 
and evolution, such as might be associated with a violent solar flare. Applications to 
calculation of neutrino fluxes in astrophysical beam dumps (with appropriately modified 
cascade propagation factors, G n ^) are also possible. 

In §4 we gave an example of the use of the semi-analytic approach in estimating "theo- 
retical" uncertainties in an expected neutrino-induced signal. A complementary approach 
is simply to compare different, independent calculations of the spectrum of atmospheric 
neutrinos. Fig. || shows a comparison of four calculations. The calculations of Ref. §| 
and [25] are one-dimensional Monte Carlo calculations. The calculation of Ref. || is an 
analytic calculation assuming a pure power-law primary spectrum of nucleons and scaling 
for the inclusive cross sections. These three calculations are within a range of ~ 10% of 
each other over the relevant range of energies. Several earlier calculations are compared 
in Fig. 7 of Ref. Q and are also within this range. The three-dimensional calculation of 
Battistoni et al. [26| is shown by the points in the figure. (Because of the complexity of 
the 3-D calculation, statistics in the near vertical bin are limited at present, especially so 
at the high-energy end p7|.) The calculation of Ref. [26] assumed the same primary spec- 



trum as used in Ref. |g] . In addition, the energy range here is well into the region for which 
three-dimensional effects can be neglected. The differences are therefore presumed to be a 
result of differences in the representation of pion and kaon production. The authors of the 
calculations shown in Fig. || are in the process of comparing hadronic interaction models 
with the goal of minimizing this source of uncertainty in the calculations of atmospheric 
muons and neutrinos. 
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Figure 8: Comparison of three calculations of the vertical flux of atmospheric neutrinos 



6 Appendix: Approximate treatment of muon decay and 
energy loss 

For muons, which lose energy and which may decay in the atmosphere before reaching 
the ground, we must first multiply by the muon survival probability and take account of 
muon energy loss between production at depth X and the ground at vertical atmospheric 
depth Xq ~ 1030g / cm 2 . This is done by multiplying Eq. [l] by the muon survival proba- 
bility before integrating over the depth, X, of pion (or kaon) decay. The muon survival 
probability given by Lipari ||] can be rewritten in the form 

P,(E,,X,X ) = ^ " ~ , „J , (14) 

\ X E fl + a(X / cos 9 - X) J 

Here E^ is the muon energy at the ground, whereas its energy at production in Eq. ^ is 
E„, = E„ + a(X /cos6 - X). 

To simplify the integration, I make the approximation of replacing X in the expression 
for E^o and in Eq. |IJ by An ~ 120 g/cm 2 , which is its mean value at low energy (E^ < 
E n <C e-x 115 GeV). Ionization energy loss is characterized by a k 2 MeV/g/cm 2 , so 
Efj,,o ~ E^ + 2 GeV/cosO, and the approximate forms will be valid for 3> 2/ cos 9 GeV, 
in practice, for E^ > 10 GeV for the applications of this paper. 



13 



For high energy, ionization energy loss is negligible, but radiative losses due to pair 
production, bremsstrahlung, and hadronic interactions of muons begin to become impor- 
tant. However, the radiation lengths are sufficiently large that the probability of significant 
radiative energy losses in the atmosphere is small, and they have been neglected here. 

With these approximations, we can now follow the same steps used to obtain Eq. |?j. 
The result is 

dN S^E^E*) 1 



Here 



and 



with 



F Nw (E w ,E N ) (15) 
(-„ \n-l 

Xn ~[ 1 + nEnCosv/en 

I A n cos a 1 \ 



dE^dE^dE^ [l — r n )E^ E n 



\ X 1 + 2/E„ cos 

C n = C»-l X ( 1 + 



(n - l)(E fl cos9 + 2)^ 
Cl =r ( 1+ (E^e + 2)) 
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